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ABSTRACT 

On large angular scales, the polarization of the CMB contains information about the evolution of 
the average ionization during the epoch of reionization. Interpretation of the polarization spectrum 
usually requires the assumption of a fixed functional form for the evolution, e.g. instantaneous reioniza- 
tion. We develop a model-independent method where a small set of principal components completely 
encapsulate the effects of reionization on the large-angle E-moAe polarization for any reionization 
history within an adjustable range in redshift. Using Markov Chain Monte Carlo methods, we apply 
this approach to both the 3-year WMAP data and simulated future data. WMAP data constrain 
two principal components of the reionization history, approximately corresponding to the total optical 
depth and the difference between the contributions to the optical depth at high and low redshifts. The 
optical depth is consistent with the constraint found in previous analyses of WMAP data that assume 
instantaneous reionization, with only slightly larger uncertainty due to the expanded set of models. 
Using the principal component approach, WMAP data also place a 95% CL upper limit of 0.08 on the 
contribution to the optical depth from redshifts z > 20. With improvements in polarization sensitivity 
and foreground modeling, approximately five of the principal components can ultimately be measured. 
Constraints on the principal components, which probe the entire reionization history, can test models 
of reionization, provide model-independent constraints on the optical depth, and detect signatures of 
high-redshift reionization. 

Subject headings: cosmic microwave background — cosmology: theory — large-scale structure of 
universe 



\ 1. INTRODUCTION 

The amplitude of the iJ-mode component of the cos- 
mic microwave background (CMB) polarization on large 
scales provides the current best constraint on the Thom- 
son scattering optical depth to reionization, r. Using 
the first three years of data from the Wilkinson Mi- 
crowave Anisotropy Probe ( WMAP) and making the sim- 
ple assum ption that the univer se was reionized instan- 

; taneously, ISpergel et all (|2006[ ) find t = 0.09 ± 0.03. 
Theoretical studies suggest that the process of reioniza- 
tion was too complex to be well described as a sudden 
transition (e.g., Barkana & Loeb 2001). Previous stud- 

I ies have examined how the constraint on r depends on 
the evolution of the globally-averaged ionized fraction 
during reionization, Xe{z), for a variety of specific theo- 
retical scenarios. If the assumed form of Xe{z) is incor- 
rect, the estimated value of r can be biased; this bias 
can be lessened by considering a wider variety of reion- 
ization histo ries at the expense of increasing the uncer- 
tainty in r (iKaBlinghat et all 120031 : iHolder et aP 120031 : 
[Colombo et al.ll2005n . 

The angular scales on which CMB polarization from 
reionization is correlated depend on the horizon size at 
the redshift of the free electrons: th e higher the red- 
shift, the higher th e multipole, i (e.g.. IZaldarriagalll997t 
iHu &: White! Il997[) . Varying Xe{z) changes the relative 
contributions to the polarization coming from different 
redshifts, and therefore changes the shape of the large- 
scale i?-mode angular power spectrum, Cf^. Because 
of this dependence, measurements of the low-^ i?-mode 
spectrum should place at least weak constraints on the 



global reionization history in addition to the constraint 
on the total optical depth. Recent studies suggest that 
WMAP data provide little information about Xe{z) be- 
yond r (e.g., Lewis ct al. 2006), but it is worth asking 
what we can ultimately expect to learn about reioniza- 
ti on from CMB polar ization. 

iHu fc Holder! ()2003l ) proposed using a principal compo- 
nent decomposition of the ionization history to quantify 
the information contained in the large-scale E'-mode po- 
larization. The effect of any ionization history on the E- 
modes can be completely described by a small number 
of eigenmode parameters, unlike a direct discretization 
of Xe(z ) in redshift bins. Here we extend the methods of 
!Hu fc H older (2003) using Markov Chain Monte Carlo 
techniques to find constraints on the principal compo- 
nents of Xe(2;) using both the 3-year WMAP observations 
and simulated future data. 

Analytic studies and simulations indicate that reion- 
ization is an inhomogencous process, and this inho- 
mogeneity is expected to contribute to the small-s cale 
CMB temperatii re an d polarization anisotropies fe.g..!Hu[ 
2000; Ihev et al.ll2dollMortonson fc Hu![2007l : iDore et all 
i2007i) . Here we focus on the large-scale i?- modes only 
{£ < 100) where such effects can be neglected, so we only 
consider the evolution of the globally-averaged ionized 
fraction as a function of redshift. 

In the following section, we describe the principal com- 
ponent method for parameterizing the ionization history 
and show that the effects of Xe{z) on large-scale £'-mode 
polarization can be encapsulated in a small set of param- 
eters. The method allows Xe{z) to be a free function of 
redshift that is not bounded by physical considerations, 
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Fig. 1. — B-mode polarization angular power spectra and ioniza- 
tion histories (inset) for a nearly-instantaneous reionization model 
with optical depth r = 0.105 [thick curves) and an extended, "dou- 
ble" reionization history with t = 0.090 (th in). Points with e rror 
bars represent the 3- year WMAP data from lPaee et al] II2006I ). 



SO in §[3] we derive limits that can be placed on the model 
parameters to eliminate most of the unphysical models 
where > 1 or aig < 0. We outline several ways to ap- 
ply the principal component approach to constrain the 
reionization history with large-scale i?-modc data in § |4l 
Using Markov Chain Monte Carlo methods, we examine 
some of these applications in more detail in § [S] using 
both 3-year WMAP data and simulated future CMB po- 
larization data. We summarize our findings and conclude 
in § ini 

2. IONIZATION HISTORY EIGENMODES 

Models with similar total optical depth but different 
ionization histories can produce markedly different pre- 
dictions for the i?-mode power spectrum. In an instanta- 
neous reionization scenario, the contribution to the op- 
tical depth from Xe{z) is concentrated at the lowest red- 
shifts possible for a given r, and the i?-mode power spec- 
trum for such a model is sharply peaked on large scales, 
at £ < 10. The main effect of shifting some portion of 
the reionization history to higher redshifts while keeping 
r fixed is to reduce the £^-mode power on the largest 
scales and increase it on smaller scales. 

This redistribution of power is illustrated in Figure [1] 
by the ionization histories and Cf^ for two models, one 
with nearly instantaneous reionization and the other with 
comparable optical depth but with Xeiz) concentrated at 
higher z. In general, a flatter large-scale E-mode spec- 
trum with power extending out to £ ~ 10-20 is a sign of 
a large ionized fraction at high redshift. However, there 
is not a one-to-one correspondence between i and z; the 
ionized fraction within any particular narrow redshift bin 
affects Cf^ over a wide range of angular scales. Ionized 
fractions in adjacent redshift bins have highly correlated 
effects on Cf^, which makes it difficult to extract from 
£^-mode polarization data a constraint on Xe at a specific 
redshift. 

As suggested bv lHu fc HoTde3 (|2003f) . one can use prin- 



cipal components of the reionization history as the model 
parameters instead of Xe in redshift bins. These com- 
ponents are defined to have uncorrelated contributions 
to the .E-mode power; since each has a unique effect 
on Cf^, the amplitudes of the components can be in- 
ferred from measurements of the large-scale power spec- 
trum. Principal component analysis of Xe{z) also indi- 
cates which components can be determined best from the 
data. In this section, we describe the principal compo- 
nent method and introduce the notation that we will use 
throughout the paper. 

Consider a binned ionization history Xe{zi), i £ 
{1, 2, . . . , iV,}, with redshift bins of width Az span- 
ning Zmin < z < Zniax, whcrc zi = z^in + Az and 

ZN, = ^max - Az so that A^, -I- 1 = (Zii,ax " Zniin)/^Z. 

Throughout this paper we assume that the ionized frac- 
tion is Xe ~ 1 for redshifts z < Zmin and Xe ~ at 



z > z„ 



We take z^ 



6, consiste nt with observa- 



tions of quasar spectra ()Fan et al.ll2006D . 

The principal components of Xe{zi) are eigenfunctions 
of the Fisher matrix, computed by taking derivatives of 
Cf^ with respect to Xe{zi): 



1=2 



dXe(zi) dXe{Zj) 



(1) 



assuming full sky coverage and neglecting noise. Since 
Xe{zi) only significantly contributes to the i?-mode spec- 
trum at small we typically truncate the sum in equa- 
tion ll]) at -^max = 100 where Cf^ is dominated by 
the first acoustic peak. The derivatives are evaluated 
at a fiducial reionization history, x^'^{zi). Following 
iHu k Holdel (|2003l ) we typically choose xl'^{zi) to be 
constant during reionization, although other functional 
forms may be used. 

Since the effects on Cf^ of Xe{zi) in adjacent redshift 
bins are highly correlated, the Fisher matrix contains 
large off-diagonal elements. The principal components 
S^(zi) are the eigenfunctions of Fij, 



F,, = (TV, + 1) 



^=1 



S,.{zi)<J^'^S^,{z,), 



(2) 



where the factor (A, -I- 1)^^ is included so that the eigen- 
functions and their amplitudes have certain convenient 
properties. The inverse eigenvalues, cr^, give the esti- 
mated variance of each principal component eigenmode 
from the measurement of low-£ i?-modes. We order the 
modes so that the best-constrained principal components 
(smallest cr^) have the lowest values of /i, starting at 
11=1. The noise level and other characteristics of an 
experiment can be included in the construction of the 
eigenfunctions, but since the effect on S^j,{z) is small we 
always use the noise-free eigenfunctions here. 

The eigenfunctions satisfy the orthogonality and com- 
pleteness relations 



dz Sf^{z)Sy{z) = (z,nax - Z^in)S^ 



(3) 



Y,S^{z,)S^{z,) = {N, + l)5,,. (4) 

The normalization of Sfj_{z) is chosen so that the eigen- 
functions are independent of bin width as Az — s- 0. In 
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Fig. 2. — Eigenfunctions for fiducial models with Zmax = 30 and 
constant Xe = 0.15 (top), and Zmax = 20 and Xe = 0.3 (bottom). 
In each case, the first three modes are shown: fi = 1 (thick), 2 
(medium) and 3 (thin). For both fiducial models, ^min = 6 and 
Az = 0.25. 

equation ([3]) and elsewhere in this paper where there are 
sums over redshift, we assume the continuous limit, re- 
placing J2i by / dz. As long as the bin width is chosen 
to be sufficiently small, the final results we obtain are in- 
dependent of the redshift binning. We adopt Az = 0.25 
as the default bin width. 

The three lowest- variance eigenfunctions for two dif- 
ferent fiducial models are shown in Figure O The lowest 
eigenmode {fi — 1) is an average of the ionized fraction 
over the entire redshift range, weighted at high z. The 
jj, = 2 mode can be thought of as a difference between 
the amount of ionization at high z and at low z, and 
higher modes follow this pattern with weighted averages 
of Xe{z) that oscillate with higher and higher frequency 
in redshift. Eigenfunctions of fiducial models with differ- 
ent values of Zmax have similar shapes with the redshift 
axis rescaled according to the width of (zmax — ^^min)- 
The eigenfunctions are mostly insensitive to the choice 
of ionized fraction in the constant-a;e fiducial histories. 

An arbitrary reionization history can be represented in 
terms of the eigenfunctions as 

Xe{z) = x^'^{z) -\-'^m^Sf^{z), (5) 
where the amplitude of eigenmode for a perturbation 

5Xe{z) = Xe{z) — X^'^{z) is 

1 /•^■°- 

= / dz S^,{z)Sxe{z). (6) 

^max ^min J Zj^jn 

[Note that our conventions f or the normahzation of to^ 
and Sf,{z) differ from those of lHu & Hold"e^ (pOO^h by fac- 
tors of {Nz + 1)^^^-] Any global ionization history Xe{z) 
over the range Zmin < z < z^ax is completely specified 
by a set of mode amplitudes m^. 

If perturbations to the fiducial history are small, 
Sxe{z) 1, then the mode amplitudes are uncorrelated, 
with covariance matrix {m^mj^) = For a fixed 

fiducial model, however, arbitrary reionization histories 
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Fig. 3. — Change in \nCf^ per unit for fiducial Xe(z) with 
Zmax = 30 (solid) and 20 (dashed), showing the effect on Cf^ of 
the first two principal components, /J. = I (thick) and 2 (thin). 

generally have Sxe ~ 1, in which case the amplitudes of 
different modes can become correlated as we discuss in 

The main advantage of using the principal component 
eigenmodes of Xe{z) instead of some other parameteriza- 
tion is that most of the information relevant for large- 
scale i?-modes is contained in the first few modes. This 
means that if one constructs Xe{z) from equation (O 
keeping only the first few terms in the sum over /i, then 
the i?-mode spectrum of the resulting ionization history 
will closely match that of Xe{z) with all modes included 
in the sum. The effect of each eigenmode on Cf^ be- 
comes smaller as fi i ncreases, as shown in Figure [3] for 
the first two modes. |Hu & HoldeaJ (|2003l ) demonstrated 
that for a specific fiducial ionization history Xg'^(z) and 
assumed true history, only the first three modes of Xe{z) 
are needed to produce Cf^ indistinguishable from the 
true i?-mode spectrum. 

The ionization histories and corresponding i?-mode 
spectra in Figure |3] demonstrate this completeness for 
a fairly extreme model in which the first ten eigenmodes 
all have significant amplitudes. Even in the simplest case 
where a single eigenmode is used in place of the original 
Xe{z), the error in Cf^ is only ~ 10%. With 3-5 modes, 
the error is a few percent or less at all multipoles and 
safely smaller than the cosmic variance of 

The top panel of Figure d] shows that this completeness 
does not extend to the ionization history itself: Xe{z) 
constructed from as many as five eigenmodes is a poor 
approximation to the full reionization history. 

From this and other similar tests on the completeness 
m Cf ^ of the lowest -variance eigenmodes we conclude 
that the first 3-5 modes contain essentially all of the 
information about the reionization history that is rele- 
vant for large-scale £^-mode polarization. This fact is 
particularly useful for constraining the global reioniza- 
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Fig. 4. — Test of completeness in for ionization histories 

with truncated sets of eigenfunctions. The base history (top panel, 
dashed) is the sum of the first 10 eigenfunctions with comparable 
values of |m(j| for each mode. Solid curves show Xe{z) (top) and 
the error in Cf^ (bottom) for histories that retain the first one 
(thin curves), three (medium), and five (thick) eigenmodes. The 
fiducial model used for the eigenmode decomposition has constant 

Xe = 0.15 from Zmin = 6 up to Zmax = 30. 

tion history with CMB polarization data using Markov 
Chain Monte Carlo techniques. Since the number of pa- 
rameters that must be added to a Monte Carlo chain 
to describe an arbitrary Xe{z) is relatively small, we can 
obtain constraints from the data that are independent 
of assumptions about the reionizatio n history with min- 
imal added computational expense ( Hu fc Hol der 2003). 
The exact number of modes required varies depending 
on the true ionization history and the fiducial history, so 
when analyzing data it is a good idea to check that the 
results do not change significantly when the next modes 
are included in the sum in equation ([5|). 

The main caveat to this model independence is that in 
practice we must set some maximum redshift for histories 
in any particular chain of Monte Carlo samples, ignoring 
any contribution to the observed low-£ Cf^ from ioniza- 
tion at z > Zniax- Since the eigenfunctions of the ion- 
ization history are stretched in redshift as Zmax increases 
(Fig. [2]), at higher Zmax more modes are needed to ac- 
curately represent any particular feature in Xe{z). For 
example, take the true ionization history to be instanta- 
neous reionization at z = 11.5. Figure [5] shows the error 
in Cf^ if we truncate the eigenmode sum of equation ([5]) 
at A'" modes using fiducial histories with z^ax = 30 and 
•Zmax — 20. For each fiducial model, the error decreases 
as the number of modes in the sum increases, but the 
error at fixed N is larger for Zmax = 30 than z,„ax = 20. 
The requirement of retaining a larger set of parame- 
ters as Zmax increases makes it less practical to study 
models with significant reio nization at extreme l y hig h 
redshifts fzmax > 1 00; e.g., iNaselskv fc Chiand 120041 ): 
iKasuva et al.l (|2004[ )]. but even for z^ax as high as ^ 40 
the number of eigenmodes needed is reasonably small 
(^<5). 

While theories of reionization provide useful priors on 



Fig. 5. — Error in Cf^ from Xe(z) constructed from the first 
N eigenmodes of fiducial models with Zmax = 30 (solid curves) 
and 20 (dashed), taking N = 3 (thin) and N = 5 (thick). The 
true Xe(z) and Cf^ are assumed to be those of an instantaneous 
reionization model with r = 0.09. 

•Zmax in the context of specific models, it would be bet- 
ter to be able to constrain z^ax empirically by mea- 
suring the i?-mode power accurately up to ^ ~ 100. 
The current 3-year WMAP polarization data have high 
enough signal-to-noise to be useful for parameter con- 
straints only a,t £ < 2A (Page et al. 2006), so their sensi- 
tivity to high-z reionization is limited. Given that ion- 
ization at some redshift generates polarization out to a 
maximum £, WMAP data can still place weak bounds 
on the total optical depth contribution above a certain 
redshift, even if it cannot distinguish whether these con- 
tributions arise from redshifts above a chosen Zmax as we 
discuss in § 15.11 The values we choose here for Zmax are 
partly influenced by the fact that few theoretical reion- 
ization scenarios predict an ionized fraction at z > 30 
that would significantly affect large-scale Cf^ . Future 
data should better constrain Cf^ at higher multipoles, 
allowing useful limits to be placed on Zmax from the data 
alone. 

The need to consider a limited range in redshift is 
shared by other methods, for example those that con- 
strain binned xAz j) instead of the eigenmode amplitudes 
(|Lewis et al.ll2006[ ). As already mentioned, the principal 
component approach has the unique advantages that re- 
sults can be made independent of bin width and that rel- 
atively few extra parameters are required. However, this 
approach also has a unique difficulty in that the physical- 
ity of the reionization history [0 < Xe{z) < 1] is not built 
in to the method. Hence constraints derived from mea- 
surements of the eigenmodes can be weaker than those 
from a method that enforces physicality. We can, how- 
ever, place some prior constraints on the set of eigenmode 
amplitudes that must be satisfied by any physical model, 
as we discuss in the next section. 

3. PRIORS FROM PHYSICALITY 

Although the actual ionized fraction must be between 
and 1 (neglecting helium reionization and the small 
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residual ionized fraction after recombination), there is 
nothing in the construction of Xe{z) from the eigenmodes 
in equation ^ that ensures that the ionized fraction will 
obey these limits. Whether or not Xe has a physical 
value at a particular redshift depends on the amplitudes 
of all of the principal components. Even for a physical 
ionization history, the truncated sum up to mode N, 



N 



„fid 



(^) 



(8) 



is not necessarily bounded by and 1 at all redshifts. 
While formally it is possible to evaluate Cf^ for reioniza- 
tion histories with unphysical values of Xe, we would like 
to eliminate as much as possible those models for which 
the full sum of the eigenmodes, Xe{z) = ImiM^oc xi^\z), 
is unphysical. 

We find the largest and smallest values of that are 
consistent with Xe{z) G [0, 1] for all z using the definition 
of in equation We are free to choose x^'^{z) e 
[0, 1] so that -x^'^{z) < Sxeiz) < 1 - x^'^iz), where the 
lower limit is strictly negative or zero and the upper limit 
is positive or zero. For a particular mode fi, the choice 
of 6xe{z) that maximizes is 

-x^d(z), < 

l-xfd(z), >0. 



Sxr'-iz) 



(9) 



Using this in equation ([6|) gives an upper limit on m^. 
Similarly, a lower limit can be obtained by reversing the 
signs of the inequalities in equation ([9]). The resulting 

physicality bounds are 

S^{z)[l 



dz- 



where 

2xl\z)]±\S^{z)\ 



2{z„ 







(10) 



If violates these bounds for any /i, the reionization 
history is guaranteed to be unphysical for some range in 
redshift. The opposite is not true, however: even if all 
satisfy equation (|10p . Xf^{z) may still be unphysical 
for some z. 

The parameter space that physical models may oc- 
cupy is restricted further by an inequality that must be 
satisfied by all eigenmodes simultaneously. Assume for 
simplicity that the fiducial model has a constant ionized 
fraction, x^*^ G [0, 1], for Zmin < z < Zmax- Any physical 
reionization history Xe{z) must satisfy 



dz[Xe{z) - X^'^Y < (Z, 



0/, 



(11) 



where / = max [(a 



.(1 



„fid\2l 



Using equation (O, 



the left side of the inequality can also be written 



dz[Xe{z) ~ x'^^'Y ^ 



dz 



mf_,Sf_,{z) 



- [Zn 



J2 K ■ 



(12) 



where the second line follows from the orthogonality of 
the eigenfunctions (eq. [3]). Comparing equations (jlip 
and (I12|) we obtain a constraint on the sum of the squares 
of the mode amplitudes, 



</, 



(13) 



where 0.25 < / < 1, depending on the value of x^'^. 
As with the physicality bounds of equation (fTO|) . this 
upper limit is a necessary but not sufficient condition for 
physicality. 

Since in practice we can only constrain a limited set 
of eigenmodes, the uncertainty in modes higher than the 
first few prevents us from simply excluding all models 
where Xg < or Xe > 1 at any redshift because the higher 
modes can have a significant effect on Xe{z). However, as 
shown in §[21 the higher modes do not affect the polariza- 
tion power spectrum since the high-frequency oscillations 
in redshift of higher modes are averaged out. Similarly, 
such eigenmodes have a small effect on the optical depth 
from a sufficiently large range in redshift. For example, 
the optical depth from 15 < z < 30 due to modes /j, > 5 
subject to the physicality constraints of this section can 
be no larger than ~ 0.01, and is likely to be smaller for 
realistic reionization scenarios. Because of this, we as- 
sume that we can place priors on the optical depth that 
correspond to < Xe < 1 over the relevant range in red- 
shift. Reionization histories with an unphysical optical 
depth over a large range in redshift are considered to be 
unphysical models since the addition of higher modes can 
not perturb the optical depth enough to give it a physical 
value. 

4. APPLICATIONS OF THE PRINCIPAL COMPONENT 
METHOD 

Once constraints on the principal components of the 
reionization history have been obtained from CMB po- 
larization data, there are several ways to use those con- 
straints to place limits on obscrvablcs such as r or to test 
theories of reionization (Hu fc Ho lder 2003). We describe 
some possible applications in this section, and in § [5] we 
put these ideas into practice using the 3-year WMAP 
data and simulated future data. 

As mentioned in § [l] the constraint on the total opti- 
cal depth to reionization depends on the assumed model 
for Xe{z). The principal component method allows us to 
explore all globally-averaged ionization histories within 
a chosen redshift range, Zmin < z < Zmax- For a given 
set of eigenmode amplitudes, {m^}, equation ([5]) yields 
the corresponding ionization history which can then be 
integrated to find the optical depth between any two red- 
shifts zi and Z2, 

r(zi,Z2) = 0.0691(1 -Fp)r!fc/i [ ' dz Xe{z). 

Jzi ^{z)/Hq 

(.14) 

In particular, the total optical depth to reionization is 



(15) 



where t(0, Zmin) ~ 0.04 for 2„iin = 6. The principal com- 
ponent approach provides a model-independent way to 
constrain t, so we expect the results to be unbiased and 
the uncertainty in r to accurately reflect the present un- 
certainty about Xe{z). Moreover, the information about 
''■(^mim ^max) IS encapsulated in the first few eigenmodes, 
i.e. the truncated ionization history of equation ([5]), 

xi^\z), with N = 3-5 for typical fiducial models. 

The total optical depth and the first principal compo- 
nent amplitude, mi, are similar quantities in that they 
are both averages of Xe{z) weighted at high z. As we 
show in the next section, CMB i?-mode polarization 
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data can constrain higher eigenmodes as well > 2). 
In particular, the second mode should be the next best 
constrained quantity since it is constructed to have the 
smallest variance after mi. Since m2 is related to the dif- 
ference in Xe between high redshift and low redshift (see 
Fig. [2]), we might guess that besides total r the data are 
mainly constraining the fraction of optical depth com- 
ing from high redshift versus low redshift. Given any 
set of mode amplitudes defining an ionization history, 
we can compute a low-z optical depth, r(zmin, ^^mid), and 
a high-z optical depth, T{zaiid, ^max), for some choice of 
intermediate redshift, Zmid- 

Note that if either of the redshift ranges [zmin, ^mid] or 
[zniid , ^^max] is too narrow, constraints on the partial op- 
tical depths will be influenced by the physicality priors 
(since < Xe < 1 sets limits on r in redshift intervals as 
discussed in § [3]) and by the uncertainty in modes higher 
than those included in the chain, which have greater ef- 
fect on the optical depth in narrower redshift intervals. 
For these reasons, Zmid should be chosen to be not too 
close to either Zmin or z,„ax- 

Given an appropriately chosen value of Zmid, con- 
straints on the principal components can be converted 
into constraints on the optical depths at high and low 
redshift. These partial optical depths are observables in 
the sense that high-z and \ow-z optical depth affect the 
large-scale E'-modes over different ranges of multipoles. 
For example, compare the two reionization models in Fig- 
ure [H Both have similar total optical depth, but in one 
case T only comes from z < 15, resulting in more power 
at £ < 8 and less power at ^ > 8 than the other model in 
which the optical depth primarily comes from z > 15. 

Besides learning about the relative amount of ioniza- 
tion at z < Zmid and z > Zmid, the partial optical depth 
constraints also provide a way to empirically set the max- 
imum redshift, Zniax- If some set of data are found to 
place a tight upper bound on r(zmid, Zmax) for fairly con- 
servative (i.e. large) values of Zmid and Zmax, then for 
analyses of future data this value of z^id can be used 
as a new, lower value for Zmax since optical depth from 
higher redshifts is small. This approach assumes that 
there is essentially no reionization earlier than the orig- 
inal Zmax, but as long as this initial maximum redshift 
is taken to be large the presence or absence of high-z 
ionization can be tested in this way over a wide range of 
redshifts. We explore this idea further in § [5j 

While we do not examine constraints on specific reion- 
ization models in this paper, the principal component ap- 
proach is well suited to model testing. Consider a model 
of the global reionization history, Xe{z] 0), parameterized 
by 6. This could be a simple toy model (for example, in- 
stantaneous reionization where is a single parameter, 
the redshift of reionization) or a more physical model 
where 6 might include parameters that govern proper- 
ties of the ionizing sources. For a particular choice of 
Xg'^(z), the principal component amplitudes of the reion- 
ization model follow from equation ([5]), giving a set of 
mode amplitudes that depend on the model parameters, 
{mp(0)}. Then constraints on the mode amplitudes from 
CMB polarization data (defined using the same fiducial 
history) can be mapped to constraints on the parameters 
of the reionization model. 

Applying this method to a model that has as one of 
its parameters a maximum redshift allows one to obtain 



constraints on Zmax within a class of theoretical models. 
Although obtaining good constraints on model parame- 
ters does not necessarily imply validation of the model 
class, this procedure provides a straightforward method 
by which different model classes can tested and falsified 
within a single analysis. 

Generating Monte Carlo chains to find constraints on 
{m^} as described in the next section can be a some- 
what time-consuming process, but it only needs to be 
done once per redshift range, after which any model of 
the global reionization history within this range Zmin < 
z < Zmax can be tested using the same parameter chains. 
In cases where constraints from the data turn out to be 
close to Gaussian, the covariance matrix of the principal 
components can be used in place of the full Monte Carlo 
chains, reducing the amount of information needed for 
model testing from ~ 10^ numbers in chains of Monte 
Carlo samples to only N{N + 3)/2 numbers when N 
eigenmodes are included in the chains [N mean val- 
ues plus N{N + l)/2 entries in the covariance matrix, 
{mfj_m^)]. However, near-Gaussian constraints on {m^} 
are likely to be possible only for certain realizations of 
future data, as discussed in the next section. 

5. MARKOV CHAIN MONTE CARLO CONSTRAINTS ON 
EIGENMODES 

We find reionization constraints from CMB polariza- 
tion data using Markov Chain Monte Carlo (MCMC) 
techniques to e xplore the principa l component parameter 
space (see e.g. IChristenseii et al.l 120011 : iKosowskv et al.l 
[200l iDunklev et all 120051) . Chains of Monte Carlo 
samples are generated using the p ublicly available code 
CosmoMC^ CLewis & Bridle^ '2002f), which includes the 
code CAMB (Lewis ct al. 2000) for computing the the- 
oretical angular power spectrum at each point in the 
{mi, . . . , mjv} parameter space. We have modified both 
codes to allow specification of an arbitrary reionization 
history calculated from a set of mode amplitudes using 
equation ([5]). 

The principal component amplitudes are the only pa- 
rameters allowed to vary in the chains. Since nearly all 
the information in Cf from reionization is contained 
within the first few eigenmodes of Xe{z), we include the 
modes {rrii, . . . , m^r} in each chain with 3 < iV < 5. 
We use only the i?-mode polarization data for parame- 
ter constraints, and assume that the values of the stan- 
dard ACDM parameters (besides r) are fixed by mea- 
surements of the CMB temperature anisotropies. This 
leads to a slight underestimate of the error on r as we 
discuss later, but to a good approximation the effect of 
Xe{z) on the large-scale £^-modes is independent of the 
other parameters. 

Specifically, we take Qbh"^ = 0.0222, flch"^ = 0.106, 
1006* = 1.04 (corresponding to h ^ 0.73), A^e"^^ = 
1.7 X 10~^, and Us = 0.96, consistent with results from 
the most recent version of the WMAP 3-year likelihood 
code. When computing the optical depth to reionization 
we take the primordial helium fraction to be = 0.24. 
We also assume that Zmin = 6, so that the optical depth 
contributed at lower redshifts is fixed at T(zmin) ~ 0.04. 
The remaining total optical depth from reionization, 
T(zmin, ^^max), IS determined by the values of {m^} for 

^ http:/ /cosmologist.info/cosmomc/ 
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each sample in the chains. The defauh bin width for 
our fiducial models x^'^{z) is Az — 0.25, which is small 
enough that numerical effects related to binning should 
be negligible. 

To get accurate results from MCMC analysis, it is im- 
portant to make sure that the parameter chains contain 
enough independent samples covering a sufficient volume 
of parameter space so that the density of the samples 
converges to the actual posterior probability distribution. 
For each scenario that we study, we run 4 separate chains 
until the Gelman and Rubin convergence statistic R, cor- 
responding to the ratio of the variance of parameters be- 
tween chains to the variance wi thin each chain, satisfies 
R-l< 0.01 (jGelman fc Rubini ri992: Brooks & Gelmaii 
[1998^. The convergence diagnostic of iRafterv fc LewisI 
(jl992) is used to determine how much each chain must 
be thinned to obtain independent samples. Both of these 
statistics and other diagnostic measures are computed 
automatically by CosmoMC. 

Since Xe{z) during reionization must match onto ~ 
1 at z = Zmin and Xe ~ at z = Zmax, there are often 
sharp transitions at these redshifts since nothing in equa- 
tion ([5|) forces Xe(z) to satisfy these boundary conditions. 
To avoid problems with the time integration in CAMB, 
we smooth the reionization history by convolving Xe{z) 
with a Gaussian of width = 0.5. 

As described in previous sections, we assume that 
< Xe < i between Zmin and z^ax- This assumption 
neglects helium reionization, which can make Xe slightly 
larger than unity, and the small residual ionized fraction 
remaining after recombination that prevents x,, from ever 
being exactly zero. These are relatively small effects, es- 
pecially since we are not placing constraints on Xe{z) di- 
rectly but rather on weighted averages of Xe over redshift. 
The residual ionized fraction at z > z^ax is accounted for 
in the Monte Carlo exploration of reionization histories. 

In § 15. H we examine the current constraints from 
the 3-year WMAP data. We then provide forecasts for 
principal component constraints with idealized, cosmic 
variance-limited, simulated data in § 15.21 In each case the 
likelihood computation includes only the £^-mode polar- 
ization data, up to ^ = 100 for simulated data and £ — 23 
for WMAP; the likehhood code for WMAP does not use 
Cf^ at smaller scales due to low signal-to-noise. At mul- 
tipoles £ > 100 the global reionization history only affects 
the amplitude of the angular power spectra, which we fix 
by setting A^e"^"^ constant in the Monte Carlo chains. A 
comparison of the MCMC results with the Fisher matrix 
approximation follows in § 15.31 

5.1. Constraints from WMAP 

In our analysis of WMAP data we use the 3-year 
WMAP likelihood code, with settings chosen so that 
likelihoods include only contributions from the low-^ 
-E-mode polarization. In this regime, the code com- 
putes model likelihoods with a pixel-based method in- 
stead of using the angular power spectrum iPage et al.l 
Since the maximum redshift at which there is still 
a significant ionized fraction is uncertain, we generate 
Monte Carlo chains using different fiducial models with 
10 < Zmax < 40. To avoid possible bias due to neglect- 
ing the possibility of high-rcdshift reionization, we focus 
here on the results obtained using the more conservative 
values of z^ax — 30 and 40. 



The MCMC constraints on the first three principal 
components from WMAP are shown in Figure [6] using 
fiducial models with Zmax = 30 and Zmax — 40. The 
marginalized constraints are plotted within the physical- 

ity bounds of equation (|10p . to|, ^ < < m^i^\ so 
the size of the contours inside each box gives an idea of 
the constraining power of the data within the space of 
potentially physical models. For both fiducial histories, 
the data place a strong upper limit on mi and weakly 
constrain m2. It is important to note that although the 
parameters {7711,7712,7713} have the same names in both 
the Zmax — 30 and Zmax = 40 plots, they are defined with 
respect to different fiducial models and so we do not ex- 
pect the contours in the left and right plots in Figure [6] 
to agree exactly. The qualitative similarity between the 
contours simply reflects the fact that the eigenfunctions 
for different z^ax have similar shapes (Figure (2). 

Since 77ii and r are both averages of Xe from Zmin to 
2,iiax with more weight at high z than low z, the strong 
constraint on r77i mostly reflects the ability of the data to 
constrain the total optical depth to reionization. Indeed, 
mi and r are strongly correlated in the Monte Carlo 
chains, while the correlations between r and higher prin- 
cipal components of Xe{z) are weaker. 

Relative to the physicality bounds, the constraints in 
Figure [6] are stronger for the Zmax = 40 chains than for 
•Zmax = 30; this is to be expected since the physicality 
bounds on {771,^} permit a greater variety of ionization 
histories, and in particular a wider range in r, when Zmax 
is larger. The range of eigenmode amplitudes allowed by 
the limits of equation (jlOp is nearly independent of Zmax, 
but the effect of a unit-amplitude principal component 
on Cf^ increases with Zmax as illustrated in Figure [3l 

As described in § 21 principal component constraints 
from MCMC yield model-independent constraints on the 
total optical depth. For each Monte Carlo sample of the 
principal components we compute r using equation (fTS]) . 
The constraints on r are listed in Table [1] for various val- 
ues of Zmax and N, the number of modes in the Monte 
Carlo chains. The uncertainty in r from the Monte 
Carlo chains is slightly underestimated because we fix 
all cosmological parameters that are not directly related 
to reionization. To estimate the effect of this assump- 
tion, we compare constraints on r in two cases where in- 
stantaneous reionization is assumed instead of using the 
model-independent principal component approach: one 
in which the non-reionization parameters are fixed as in 
the rest of our Monte Carlo ch ains, and the other, f rom 
the 3-year WMAP analysis of iSpergel erall (|2006D . in 
which these parameters are allowed to vary. These con- 
straints, listed in rows 4 and 5 of Table [U suggest that 
fixing ACDM parameters besides r reduces ar by ^ 10%. 

In all cases where we fix the other cosmological param- 
eters, the uncertainty in r is similar regardless of whether 
we use principal components to explore a variety of reion- 
ization histories or restrict the analysis to instantaneous 
Xe{z). This is somewhat surprising, since in general one 
would expect that expanding the model space for Xeiz) 
would increase the estimated error on r. The physicality 
priors may be responsible for reducing ar slightly — top- 
hat priors on {m^} induce a prior on r that is flat over a 
certain range but falls off approximately linearly at the 
edges — but even after accounting for priors, the optical 
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Fig. 6. — Marginalized 2D contours {dark shading: 68% CL; light: 95%) for the first three principal components of ^^(z) from MCMC 
analysis of 3-year WMAP E-raode polarization data, using fiducial reionization histories with Zmax = 30 (left) and Zmax = 40 (right). Note 
that the mode amplitudes are defined diff'erently in the left and right panels since the fiducial models are different. 
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Constraints on total optical depth to reionization. 
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depth constraint is robust to replacing the instantaneous 
reionization assumption with a model-independent anal- 
ysis. 

The insensitivity of the constraint on r to the set of 
models considered is partly due to the fact that the de- 
generacy in the eigenmode constraints is aligned in the 
direction of constant t, as shown in Figure [7| in the 
mi — m2 plane. The set of instantaneous reionization 
models, plotted as a curve in Figure [71 cuts across the 
WMAP constraints on more general models in a region 
of high posterior probability where the distribution of 
samples varies slowly along lines of constant r. This 
large overlap between the general reionization histories 
favored by the data and the line of instantaneous mod- 
els is the main reason why the probability distributions 
P(t), plotted in the right panel of Figure [71 are simi- 
lar for the two classes of models. [The sharp cutoff at 
low T in the instantaneous reionization P(t) comes from 
our 2„iin — 6 prior.] Note that the fact that the instan- 
taneous reionization curve passes through the middle of 
the 68% confidence region also indicates that models with 
rapid reionization are not at all disfavored by the 3-year 
WMAP data. We use two methods to find the posterior 
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Fig. 7. — Left panel: Marginalized 2D contours (68 and 95%) 
from a 2max = 30 WMAP chain of Monte Carlo samples (thick 
curves), lines of constant t = {0.06,0.09,0.12} (dashed lines), and 
the set of one-parameter instantaneous reionization models (thin 
curve) plotted in the mi — m2 plane. Higher eigenmodes (p, > 3), 
which only weakly influence r, are fixed to = for the purposes 
of plotting the constant-r lines. Right panel: posterior probability 
of the optical depth from 3-year WMAP data, P(t), for arbitrary 



Xe(z) with 



■ 6, z„ 



■ 30, and N = 5 principal components 



(thick curve), and for the one-parameter instantaneous reionization 
model with a Zmin = 6 prior (thin curve). 



distribution for r in the instantaneous reionization case: 
one is the usual approach of varying the optical depth 
(or reionization redshift) in a Monte Carlo chain, and 
the other involves computing P(t) for a subset of sam- 
ples from chains in which principal components of Xe{z) 
are varied, selecting only those samples with {m^} values 
close to the ID instantaneous reionization curve. Both 
approaches produce consistent probability distributions; 
the former method is used for P(t) plotted in Figure [71 
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Fig. 8. — Constraints on optical depth during reionization di- 
vided between two wide ranges in redshift, r(6, 20) and t(20, ^max), 
for Zmax = 30 (ie/t) and 40 (right). The contours are drawn at 68% 
{dark gray) and 95% (light gray) 2D marginaUzed posterior proba- 
biUty. Dashed lines in each panel are lines of constant total optical 
depth: r = 0.12 for the upper line, and r = 0.06 for the lower line. 

The weak TO2 constraint visible in Figure [S] suggests 
that in addition to determining the total optical depth, 
WMAP data may provide useful limits on high-z and 
low-z optical depth as defined in the previous section. 
The 68 and 95% posterior probability contours for the 
partial optical depths from WMAP are shown in Fig- 
ure [8] for Zmid = 20 and z^ax = {30,40}. The Xe ^ 
physicality prior cuts off the contours at t(zi, 22) — 0; in 
both panels, the upper limits set by a;e = 1 are outside 
the plotted area, so the contours are not strongly influ- 
enced by those priors. The error "ellipses" are narrowest 
in the direction along which the total optical depth is 
constrained, as shown by dashed lines of constant total 
optical depth at r = 0.06 and r = 0.12 [approximately 
the upper and lower 1 a hmits from the 3-year WMAP 
analysis of ISpergel et~al] (|2006l )]. 

Comparison of the two panels in Figure [5] reveals that 
the choice of z^ax does not have a large effect on con- 
straints on the optical depth above and below z = 20. 
This suggests that r(20, Zmax) should be interpreted as 
t(z > 20). We provide further justification for this in- 
terpretation in § 15.21 

The 95% upper limits on the high-redshift optical 
depth, t{z > 20), are 0.076 and 0.078 for z,„ax — 30 and 
40, respectively, after marginalizing over all other param- 
eters. This is not a particularly strong constraint, since 
this result is only marginally inconsistent with all of the 
optical depth from reionization coming from z > 20, but 
with future data it should be possible to either reduce 
the upper limit on high-z optical depth or to detect the 
presence of a substantial ionized fraction at high redshift 
(see SlOl. 

Models of reionization can be tested by computing 
the principal components of proposed ionization histories 
and comparing with constraints on {m^} from Monte 
Carlo chains. The WMAP constraints in Figure [5] are 
non-Gaussian, in part because the physicality priors on 



{m^} intersect the posterior probability contours where 
the likelihood is large. Because the constraints do not 
have a simple Gaussian form, it appears that the full pa- 
rameter chains are necessary for accurate model testing, 
although the viability of various models can be estimated 
by comparing their eigenmode amplitudes with marginal- 
ized constraints such as those in Figure [HI keeping in 
mind that models favored by marginalized constraints 
could be disfavored in the full A'^-D parameter space. 

5.2. Cosmic variance-limited data 

To forecast how well principal components of the reion- 
ization history could be measured by future CMB polar- 
ization experiments, we repeat the analysis of ii l5.1l using 
simulated realizations of the full-sky, noiseless £^-mode 
angular power spectrum instead of WMAP data. Any 
real experiment will of course involve sky cuts, noise, 
foregrounds, and other complications, so the results pre- 
sented here represent an optimistic limit on what we can 
learn about the global reionization history from low-i 
E-mode polarization. At the end of this section, we esti- 
mate how much the constraints might be degraded from 
this idealized case for an experiment with characteristics 
similar to those proposed for the Planck satellite. 

We generate simulated realizations of Cf^ drawn from 
distributions with cosmic variance determined by the 
theoretical angular power spectra that we compute using 
CAME. For the jth sample in a chain, the likelihood 
including only £'-mode polarization is 

/ 1 \ / r<EE fiEE \ 

where Cf^ is the spectrum of the simulated data and 
is the theoretical spectrum calculated with the pa- 
rameter values at step j in the chain. 

We have run Monte Carlo chains for multiple realiza- 
tions of Cf ^ drawn from spectra computed assuming a 
variety of "true" reionization histories, x^^^''{z). We start 
by taking the true history to be a model with nearly in- 
stantaneous reionization and r = 0.105. As a contrasting 
model for comparison we also use an extended, double 
reionization model with r — 0.090. Figure [T] shows Cf^ 
and Xf.{z) for each of these models. 

Since parameter constraints derived using a single draw 
of Cf^ from the underlying power spectrum may contain 
features unique to that realization, we generated Monte 
Carlo chains for 10 random realizations of the instanta- 
neous reionization power spectrum. Two of these real- 
izations are plotted as points in Figure [51 along with the 
theoretical spectrum with cosmic variance bands. The 
thick curves in Figure [H] are spectra of the best-fit mod- 
els from the Monte Carlo chains for each realization. The 
overall best-fit models are similar for the two realizations 
and both agree closely with the theoretical Cf^. The 
dashed curve in the right panel of Figure [9] shows an 
"alternative" best-fit model that we discuss later. 

The 2D marginalized MCMC constraints on principal 
components for these two realizations are shown in Fig- 
ure [TOl The plotted regions are bounded by the physical 
top-hat priors on {m^} as in Figure O but the eigen- 
modes are different here with z^ax — 20 instead of 30 
or 40 for the fiducial history. Since we have constructed 
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Fig. 9. — _E-mode spectra of the best-fit models from two Monte Carlo chains (thick curves), each using a different realization of simulated 
cosmic variance-limited data (points). The realizations are drawn from the same of the instantaneous reionization model in Figure [Tl 

plotted as thin curves with shaded regions showing the 1 cr cosmic v aria nce limits. The dashed curve in the right panel shows the best fit 
model within the high-m2 peak in the posterior probability (see Fig. IIOI I . 
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Fig. 10. — Marginalized 2D contours (68 and 95%) for the first three principal components with two different realizations of cosmic 
variance-limited data simulated using the instantaneous reionization history. The fiducial model has maximum redshift Zmax = 20, and 
only the three modes shown here are varied in the parameter chains. The cross in each panel marks the location of the principal components 
of a;t"'"<=(2). 



a;*''"°(z) to have no ionization at z > 15, we know that 
•Zmax = 20 should be a large enough maximum redshift; 
this choice of Zmax also reflects the fact that improved 
empirical determination of Zmax will allow it to be set 
at lower redshift for future data analysis, assuming that 
high-z reionization is not detected. 

The contours in Figure [TOl show that cosmic variance- 
limited data can constrain models within the physically 
allowed parameter space much better than current data 
can. The chains used for the results shown here vary the 
first 3 eigenmodes; by running chains with A'^ > 3 we 
find that cosmic variance-limited data can provide 2 a 
constraints that are tighter than the physicality bounds 



for the first 4 eigenmodes of a Zmax = 20 fiducial model 
and for at least 5 eigenmodes when Zmax = 30. 

Although the only difference in the MCMC setup for 
the left and right panels of Figure \W\ is the realization 
of Cf^ (Fig. [9|), there are some significant differences 
in the resulting eigenmode constraints, particularly in 
the TO2 — plane. This kind of variation is typical 
among the realizations of simulated data, and is related 
to whether features in the theoretical Cf^ remain in- 
tact with the inclusion of cosmic variance. The relevant 
features in the spectrum are the small oscillations in the 
low-power "valley" in Cf^ on scales 10 < ^ < 30. In the 
example shown here, the feature that matters most is the 
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bump in the theoretical spectrum at ^ w 13. In the left 
panel of Figure[9l the Cf^ realization retains this bump, 
while in the right panel the bump is washed out by cosmic 
variance. The result is that the left realization is better 
able to pick out those reionization models in Monte Carlo 
chains that reproduce the true Cf^ , leading to the tight 
constraints on the left side of Figure [TOl The realization 
in the right panel of Figure[9]lacks this important "finger- 
print" for identifying models that match the theoretical 
Cf^, so the data allow a wider variety of models as re- 
flected in the constraints on the right side of Figure [TOl 
In this case, there are two best-fit models corresponding 
to the two 68% contours separated in the m2 direction. 
The spectrum of the overall best-fit model (with smaller 
7712), which is close to the theoretical spectrum, is plot- 
ted as a thick solid curve in the right panel of Figure [51 
while the other best-fit model (at larger m2) is plotted as 
a dashed curve. These two models have significantly dif- 
ferent spectra, but because of the scatter in the random 
draw of Cf^ they are each able to fit the data better 
at some multipoles and worse at others in such a way 
that the high-m2 best-fit model is only a slightly worse 
fit than the low-m2 best-fit model. 

Among the realizations of simulated data that we ex- 
plored with MCMC methods, constraints like those in 
the right panel of Figure [10] are a fairly extreme case, 
occurring about 10-20% of the time. In general the real- 
izations form a sort of continuum between the two pre- 
sented here, with about half showing some displacement 
of the 95% contour towards larger m2 but not as much 
as in the second realization that we show as an example. 

The extent to which differences between realizations 
show up in principal component constraints depends on 
a;*''"°(z). For example, the Cf^ for the double reioniza- 
tion model (Fig.[T]) have more pronounced oscillations at 
i ~ 10-30 than the instantaneous reionization spectrum, 
so they are not as easily erased by cosmic variance or 
noise and the resulting principal component constraints 
tend to be more consistent from one realization to the 
next and more like the left panel of Figure [TOl 

The differences between constraints from various re- 
alizations of simulated data suggest the interesting pos- 
sibility that our ability to learn about the reionization 
history from large scale E'-mode polarization may ulti- 
mately depend on the luck of the draw of Cf^ at our 
particular vantage point. In an unlucky draw, cosmic 
variance can distort subtle features in the power spec- 
trum in a way that would limit constraints on reioniza- 
tion eigenmodes to be worse than we might expect from 
the Fisher approximation or a more typical draw. The 
good news is that even a realization like the one in the 
right panels of Figures 01 and [TUl would permit constraints 
that are far better than what is currently possible. It is 
also important to note that although such constraints 
are weaker than in the best-case scenario, they are still 
consistent with the true parameter values and in general 
would not lead us to rule out the true reionization history 
based on a principal component analysis of the data. 

As with the principal components of Xe{z), MCMC 
analysis indicates that constraints on the total optical 
depth to reionization would be greatly improved with 
cosmic variance-limited data: ctt ~ 0.005, down a factor 
of six from the 3-year WMAP value, cr^ « 0.03. The 
constraints on r using the realization of instantaneous 



reionization Cf^ in the left panel of Figure [H are listed 
in Table [H for two fiducial models with Zmax = 20 and 
30. Optical depth constraints from the other draw of 
Cf^ are similar, so estimates of r and ar appear not 
to be biased for atypical realizations. It makes sense 
that differences in realizations do not significantly affect 
r since the optical depth constraint comes mainly from 
the peak in Cf^ on the largest scales whereas differences 
in principal component constraints between realizations 
arise from the low-power part of the spectrum at higher 

e. 

The model-independent constraints on r are neither 
biased nor significantly weaker relative to the constraint 
obtained assuming instantaneous reionization (bottom 
row in Table [J), as is true for the WMAP data. This 
should not be surprising since we have assumed that 
^truo^^-j g^j-^ instantaneous reionization model. How- 
ever, if x^^^^lz) is in fact very different from the instan- 
taneous model, the estimate of r obtained under the as- 
sumption of instantaneous reionization will be incorrect. 
For example, if we take xl'^^°{z) to be the double reion- 
ization model of Figure [l] the model-independent ap- 
proach yields an optical depth constraint consistent with 
the true value of r = 0.090. On the other hand, MCMC 
analysis restricted to instantaneous reionization histories 
gives a significantly biased estimate, r = 0.135 ± 0.005. 
In this case, the ID curve of instantaneous reionization 
models in the space of eigenmodes never lies near the 
principal component values favored by the data, so any 
analysis that only considers such models would never find 
a good fit to the data. This illustrates the importance 
of a model-independent approach: although results from 
current data may not be significantly affected by the as- 
sumption of a specific form of Xe{z), future data will 
be sensitive enough to the reionization history that the 
choice of a specific model can greatly impact the results 
(iHolder et al]l2003l : IColombo et al.1 12005'1. 

Though constraints on the contributions to r from high 
and low z are currently limited to weak upper bounds, 
future polarization data will likely enable more defini- 
tive detections of high-z reionization if it is present, or 
tighter upper limits if absent. The left and right panels 
of Figure [TTl show these constraints from Monte Carlo 
chains using the Cf^ plotted in the left and right panels 
of Figure [51 respectively. The true ionization history has 
zero ionized fraction at z > 15, and the MCMC analysis 
places a 95% CL upper limit on t{z > 15) of ^ 0.03. 
Figure [T2l shows MCMC constraints on the same partial 
optical depths using simulated data where the true model 
is instead taken to be the extended, double reionization 
history plotted as a thin curve in the inset of Figure [1] 
In this case, most of T(z,„in, z^ax) comes from z > 15, 
and using the principal components of Xe{z) the high-z 
optical depth can be detected at a high level of confi- 
dence and measured to an accuracy of ~ 0.01. Note that 
the constraints in both Figures [11] and [121 are consistent 
with the WMAP constraints in Figure [H 

As with WMAP data, most of the constraints from 
simulated data in Figures [TTl and [T^l do not appear to 
depend strongly on the choice of z,„ax- For the instan- 
taneous reionization model used for Figure 111! this is to 
be expected since we know that the simulated Cf^ have 
no contribution from z > 20. Each set of contours in 
Figure HH is consistent with the values of t(6, 15) and 
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Fig. 11. — Same as Figure |8]for simulated cosmic variance-limited data using the two realizations of instantaneous reionization Cf^ 
shown in Figure |9] Crosses indicate the partial optical depths for the true reionization history. The redshifts chosen for z^^^ and Zmax are 
also lower here than for the WMAP analysis, reflecting the expectation that future data may permit better empirical constraints on Zmax- 



0.1 




0.02 0.04 0.06 0.02 0.04 0.06 
t(15, 20) t(15, 30) 

Fig. 12. — Same as Figure [TT] but assuming that is 
a double reionization history with Zmax = 23 (thin curves in Fig- 
ure[T]l. Crosses mark the true partial optical depths for this model; 
note that unlike in Figure lTTI t(15, 20) ^ t(15, 30) for the expected 
values here since x^J^°{z) has a nonzero ionized fraction between 
redshifts 20 and 30. 



t(15, Zmax) for The constraints in the second 

panel in Figure [TT] extend to larger high-z optical depth 
than those in the first panel because the effect of the 
physicality priors is not as strong at larger Zmax, as noted 
in §[531 the leftmost Zmax = 20 contours would stretch to 
larger t(15, 20) and smaller t(6, 15) if we did not apply 
the physicality bounds. The two panels on the right side 
of Figure [TT] show weaker constraints on the partial opti- 
cal depths because the principal component constraints 
for that realization of Cf^ are weaker. 

In the extended, double reionization model used for 
Figure [T^ reionization actually starts aX z — 23, so the 
choice of Zmax = 20 is not appropriate for this model. 



The result of such an error is that the constraint on 
t(15, 20) is inconsistent with the true value, marked by 
a cross in the left panel of Figure [T^l The 2;max = 30 
constraints in the right panel are consistent with the 
expected value, although not in perfect agreement due 
to cosmic variance. It is interesting that although the 
value chosen for z^ax hi the left panel is too low for this 
model, the constraints on t{z > 15) for Zmax = 20 and 
•Zmax = 30 are still consistent with each other. This sup- 
ports the interpretation of constraints on T(z,„id, Zmax) 
as t{z > Zinid) in §[0 

While the cosmic variance-limited simulated data are 
useful for determining how well the reionization history 
could possibly be constrained by large-scale E-mode po- 
larization measurements, it is also interesting to ask how 
well we can do with future experiments that fall some- 
what short of the idealized case that we have considered 
so far. In particular, the upcoming Planck satellite is 
expected to improve our kno wledge of the large-scale E- 
mode spectrum substantially ()The Planck CoUaborationl 
120061 ): what does this imply for constraints on Xe{z)7 
To estimate what might be possible with Planck data, 
we assume that after subtracting foregrounds a single 
foreground-free frequency channel remains for constrain- 
ing the \ow-£ i5-mode polarization. We take this to be 
the 143 GHz channel with a white noise power level of 

— 1/2 

Wp = 81 /LiK' and beam size 6'fwhm = 7.1', and we 
assume that the sky coverag e is /skv = 0-8 after cut- 



ting out the Galactic plan e (|The Planck CoUaborationl 
I2006t lAlbrecht et al.l [20061 ). We compute the likelihood 

of Monte Carlo samples using the routines provided in 
CosmoMC and analyze parameter chains with the prin- 
cipal component method as described for WMAP and 
cosmic variance-limited data. 

For various choices of x^'^{z) and x'^J^°{z) we find 
that going from full-sky cosmic variance-limited data to 
Planck-like data increases the uncertainty in Xe{z) prin- 
cipal components, tr^, and in the optical depth, CTt-, by 
roughly a factor of two. Based on these results, it should 
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Fig. 13. — Monte Car lo constraints in the mi-m2 plane from the 
left panel of Figure 1101 (68 and 95% contours) compared with the 
estimated constraints from the Fisher matrix {shaded contours). 

be possible to constrain about three eigenmodes at ^ 2 a 
within the space of physical models, and we expect r 
to be determined to an accuracy of about ±0.01. This 
is consistent with previous estimates of the Planck op- 
tical depth uncertainty when considering more general 
reionization models than instantaneous reionization (e.g., 
iHolder et al.l[2003HHu fc Holder! [2003h . 

5.3. Fisher matrix predictions versus MCMC results 

The Fisher matrix analysis of § [5] predicts that the 
principal components of Xc{z) should be uncorrelated, 
with errors given by from equation Q (in the ide- 
alized limit of full sky, noiseless observations). These 
characteristics rely on the assumption that the differ- 
ence between the true and fiducial reionization histories 
is small: fe^™<=(.z) = x*™'=(z) - xl'^{z) < 1. Clearly this 
assumption will not hold for any single fiducial history if 
one wants to consider a variety of possible true histories, 
and the consequence is that {m^rai,) differs from a^S/^i^^. 

As an example of the difference between the Fisher 
approximation and the full Monte Carlo analysis, in Fig- 
ure [13] the marginalized constraints in the mi — m2 
plane from the left panel of Figure [TO] are compared 
with the Fisher matrix error ellipses, centered on the 
true (TOi,m2) for the same choices of a:g™°(z) (instanta- 
neous reionization with r = 0.105) and x^'^{z) (constant 
Xe = 0.15 out to Zmax = 20). The MCMC constraints 
have correlated values of mi and m2 with somewhat 
larger uncertainties than the Fisher matrix {cr^}. The 
eigenfunctions, Sfj_{z), are constructed to have orthogo- 
nal effects on Cf^ in the vicinity of the fiducial model, 
but for large Sxe{z) the orthogonality breaks down. For- 
tunately the errors on r remain largely unaffected be- 
cause the mi — m2 correlation is oriented such that the 
increased errors are along the direction of constant total 
r (see Figure [7]). 

One can compute new eigenfunctions by evaluating 
dlnCf^ / dxe{z) at x'^^^'^{z) and see that their depen- 
dence on redshift differs from the original Sf^{z). It is 



possible to decorrelate the principal components by di- 
agonalizing the covariance matrix from the Monte Carlo 
chains and rotating the eigenmodes into the new basis. 
However, this procedure is unnecessary for most appli- 
cations; the most important property of the principal 
components for constraining Xe{z) is that the first few 
modes form a complete basis for Cf^ on large scales, 
accurate within cosmic variance limits. Completeness of 
the eigenmodes holds true even if they are not exactly 
orthogonal. 

Finally, note that if the true ionization history and the 
draw of Cf^ permit constraints on principal components 
similar to those in the left panel of Figure [TO] then these 
constraints may be close enough to Gaussian so that the 
covariance matrix would be sufhcient to describe the in- 
formation in the Cf^ for finding constraints on reion- 
ization models. However, if what we see looks more like 
the right panel of Figure [TO] then the full chains of Monte 
Carlo samples may be necessary for further applications 
even if we can measure the i?-mode polarization per- 
fectly. 

6. DISCUSSION 

Observations of the large-scale iJ-mode polarization of 
the CMB in the near future are expected to yield new in- 
formation about the spatially-averaged reionization his- 
tory of the universe. The principal components of the 
reionization history are a promising tool for extracting 
as much of that information as possible from the data. 
We have shown that the principal component method 
can be usefully applied to real, currently available data, 
and forecasts from simulated data suggest that there is 
room to substantially improve constraints on the reion- 
ization history using this method as measurements of the 
large-scale E'-modes improve. 

We find that the key featur es of the princ i pal co mpo- 
nent analysis put forward by iHu fc Holder! (|2003l ) con- 
tinue to apply when we go from the Fisher matrix ap- 
proximation to an exploration of the full likelihood sur- 
face using Markov Chain Monte Carlo methods. For 
fairly conservative choices of the maximum redshift of 
reionization (zmax ~ 30-40), only the first five principal 
components at most are needed for a complete represen- 
tation of the i?-mode angular power spectrum to within 
cosmic variance. To account for arbitrary reionization 
histories in the analysis of CMB data, only a few addi- 
tional parameters must be included in chains of Monte 
Carlo samples if those parameters are taken to be the 
lowest- variance principal components of Xe{z). 

Specific models of reionization can be tested easily 
by computing their eigenmode amplitudes and compar- 
ing with constraints on the eigenmodes from the data. 
Constraints on derived parameters, such as total opti- 
cal depth or the optical depth from a certain range in 
redshift, represent other applications of the MCMC con- 
straints on principal components of Xe{z). 

Often, estimates of the optical depth to reionization are 
computed assuming instantaneous reionization or some 
other simple form for Xe{z). Here we extend the analysis 
of the 3-year WMAP data to allow a more general set of 
models of the global reionization history. We find that 
expanding the model space does not significantly widen 
the uncertainty in r beyond the instantaneous reioniza- 
tion value of Ut — 0.03. Robust t constraints are impor- 
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tant for tests of the dark energy based on the growth of 
structure since they control the uncertainty on the am- 
phtude of the initial spectrum. 

Moreover, even with current data the principal compo- 
nent constraints are beginning to show the possibility of 
determining properties of reionization in addition to r. 
By comparing the optical depth from high z with that 
from low z, we obtain an upper limit on the contribution 
to the optical depth from high redshift: r(z > 20) < 0.08 
at 95% confidence, assuming that there is no significant 
episode of reionization at z 3> 40. 

Due to the limitations of noise and foreground con- 
tamination, only the first two eigenmodes of the reioniza- 
tion history, toi and m2, can be determined with present 
polarization data to any reasonable degree of accuracy. 
Constraints on these two modes come primarily from the 
main, broad peak in i?-mode power at low i from reion- 
ization. 

As measurements of the i?-mode polarization improve, 
for example from additional WMAP data or through 
planned future experiments such as Planck, better knowl- 
edge of the low-power "trough" in Cf^ between the main 
reionization peak and the first acoustic peak should en- 
able constraints on the third and higher principal com- 
ponents, up to about 7715 for near cosmic variance-limited 



data. Since constraints on these higher modes rely on the 
ability to identify subtle features in the trough of Cf^, 
the ultimate accuracy to which the eigenmodes can be 
determined may depend on whether or not the necessary 
features are well reproduced in the particular random 
draw of Cf^ that is available to us. However, even if we 
are unlucky enough to have a realization in which some 
of the important features of the spectrum are washed out 
by randomness, it should still be possible to measure sev- 
eral of the principal components to better accuracy than 
is currently possible. Knowledge of the /i > 3 eigenmodes 
of Xe{z), along with improved constraints on mi and 1712, 
will allow more stringent tests of reionization models and 
a better understanding of the global reionization history. 
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